{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **ATLAS/ICESat-2 Land Ice Height [ATL06](https://nsidc.org/data/atl06/) Exploratory Data Analysis**\n",
    "\n",
    "[Yet another](https://xkcd.com/927) take on playing with ICESat-2's Land Ice Height ATL06 data,\n",
    "specfically with a focus on analyzing ice elevation changes over Antarctica.\n",
    "Specifically, this jupyter notebook will cover:\n",
    "\n",
    "- Downloading datasets from the web via [intake](https://intake.readthedocs.io)\n",
    "- Performing [Exploratory Data Analysis](https://en.wikipedia.org/wiki/Exploratory_data_analysis)\n",
    "  using the [PyData](https://pydata.org) stack (e.g. [xarray](http://xarray.pydata.org), [dask](https://dask.org))\n",
    "- Plotting figures using [Hvplot](https://hvplot.holoviz.org) and [PyGMT](https://www.pygmt.org)\n",
    "\n",
    "This is in contrast with the [icepyx](https://github.com/icesat2py/icepyx) package\n",
    "and 'official' 2019/2020 [ICESat-2 Hackweek tutorials](https://github.com/ICESAT-2HackWeek/ICESat2_hackweek_tutorials) (which are also awesome!)\n",
    "that tends to use a slightly different approach (e.g. handcoded download scripts, [h5py](http://www.h5py.org) for data reading, etc).\n",
    "The core concept here is to run things in a more intuitive and scalable (parallelizable) manner on a continent scale (rather than just a specific region)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import glob\n",
    "import json\n",
    "import logging\n",
    "import netrc\n",
    "import os\n",
    "\n",
    "import cartopy\n",
    "import dask\n",
    "import dask.distributed\n",
    "import hvplot.dask\n",
    "import hvplot.pandas\n",
    "import hvplot.xarray\n",
    "import intake\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pyproj\n",
    "import requests\n",
    "import tqdm\n",
    "import xarray as xr\n",
    "\n",
    "import deepicedrain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style=\"border: 2px solid white;\">\n",
       "<tr>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3 style=\"text-align: left;\">Client</h3>\n",
       "<ul style=\"text-align: left; list-style: none; margin: 0; padding: 0;\">\n",
       "  <li><b>Scheduler: </b>tcp://127.0.0.1:43233</li>\n",
       "  <li><b>Dashboard: </b><a href='http://127.0.0.1:8787/status' target='_blank'>http://127.0.0.1:8787/status</a></li>\n",
       "</ul>\n",
       "</td>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3 style=\"text-align: left;\">Cluster</h3>\n",
       "<ul style=\"text-align: left; list-style:none; margin: 0; padding: 0;\">\n",
       "  <li><b>Workers: </b>10</li>\n",
       "  <li><b>Cores: </b>10</li>\n",
       "  <li><b>Memory: </b>201.22 GB</li>\n",
       "</ul>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Client: 'tcp://127.0.0.1:43233' processes=10 threads=10, memory=201.22 GB>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Limit compute to 10 cores for download part using intake\n",
    "# Can possibly go up to 10 because there are 10 DPs?\n",
    "# See https://n5eil02u.ecs.nsidc.org/opendap/hyrax/catalog.xml\n",
    "client = dask.distributed.Client(n_workers=10, threads_per_worker=1)\n",
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick view\n",
    "\n",
    "Use our [intake catalog](https://intake.readthedocs.io/en/latest/catalog.html) to get some sample ATL06 data\n",
    "(while making sure we have our Earthdata credentials set up properly),\n",
    "and view it using [xarray](https://xarray.pydata.org) and [hvplot](https://hvplot.pyviz.org)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the local intake data catalog file containing ICESat-2 stuff\n",
    "catalog = intake.open_catalog(\"deepicedrain/atlas_catalog.yaml\")\n",
    "# or if the deepicedrain python package is installed, you can use either of the below:\n",
    "# catalog = deepicedrain.catalog\n",
    "# catalog = intake.cat.atlas_cat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:                (delta_time: 1558437)\n",
      "Coordinates:\n",
      "    version                <U3 '003'\n",
      "    cyclenumber            <U2 '08'\n",
      "    orbitalsegment         <U2 '11'\n",
      "    revision               <U2 '01'\n",
      "  * delta_time             (delta_time) datetime64[ns] 2020-07-16T00:33:22.06...\n",
      "    latitude               (delta_time) float64 dask.array<chunksize=(50000,), meta=np.ndarray>\n",
      "    longitude              (delta_time) float64 dask.array<chunksize=(50000,), meta=np.ndarray>\n",
      "    datetime               (delta_time) datetime64[ns] 2020-07-16T00:33:22 .....\n",
      "    referencegroundtrack   (delta_time) <U4 '0315' '0315' ... '0329' '0329'\n",
      "Data variables:\n",
      "    atl06_quality_summary  (delta_time) int8 dask.array<chunksize=(50000,), meta=np.ndarray>\n",
      "    h_li                   (delta_time) float32 dask.array<chunksize=(50000,), meta=np.ndarray>\n",
      "    h_li_sigma             (delta_time) float32 dask.array<chunksize=(50000,), meta=np.ndarray>\n",
      "    segment_id             (delta_time) float64 dask.array<chunksize=(50000,), meta=np.ndarray>\n",
      "    sigma_geo_h            (delta_time) float32 dask.array<chunksize=(50000,), meta=np.ndarray>\n",
      "Attributes:\n",
      "    Description:  The land_ice_height group contains the primary set of deriv...\n",
      "    data_rate:    Data within this group are sparse.  Data values are provide...\n"
     ]
    }
   ],
   "source": [
    "try:\n",
    "    netrc.netrc()\n",
    "except FileNotFoundError as error_msg:\n",
    "    print(\n",
    "        f\"{error_msg}, please follow instructions to create one at \"\n",
    "        \"https://nsidc.org/support/faq/what-options-are-available-bulk-downloading-data-https-earthdata-login-enabled \"\n",
    "        'basically using `echo \"machine urs.earthdata.nasa.gov login <uid> password <password>\" >> ~/.netrc`'\n",
    "    )\n",
    "    raise\n",
    "\n",
    "# data download will depend on having a .netrc file in home folder\n",
    "dataset = catalog.icesat2atl06.to_dask().unify_chunks()\n",
    "print(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dataset.hvplot.points(\n",
    "#     x=\"longitude\",\n",
    "#     y=\"latitude\",\n",
    "#     c=\"h_li\",\n",
    "#     cmap=\"Blues\",\n",
    "#     rasterize=True,\n",
    "#     hover=True,\n",
    "#     width=800,\n",
    "#     height=500,\n",
    "#     geo=True,\n",
    "#     coastline=True,\n",
    "#     crs=cartopy.crs.PlateCarree(),\n",
    "#     projection=cartopy.crs.Stereographic(central_latitude=-71),\n",
    "# )\n",
    "catalog.icesat2atl06.hvplot.quickview()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data intake\n",
    "\n",
    "Pulling in all of the raw ATL06 data (HDF5 format) from the NSIDC servers via an intake catalog file.\n",
    "Note that this will involve 100s if not 1000s of GBs of data, so make sure there's enough storage!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download all ICESAT2 ATLAS hdf files from start to end date\n",
    "dates1 = pd.date_range(start=\"2018.10.14\", end=\"2018.12.08\")  # 1st batch\n",
    "dates2 = pd.date_range(start=\"2018.12.10\", end=\"2019.06.26\")  # 2nd batch\n",
    "dates3 = pd.date_range(start=\"2019.07.26\", end=\"2020.05.13\")  # 3rd batch\n",
    "dates = dates1.append(other=dates2).append(other=dates3)\n",
    "# dates = pd.date_range(start=\"2020.05.14\", end=\"2020.07.16\")  # custom batch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Submit download jobs to Client\n",
    "futures = []\n",
    "for date in dates:\n",
    "    revision = 2 if date in pd.date_range(start=\"2020.04.22\", end=\"2020.05.04\") else 1\n",
    "    source = catalog.icesat2atlasdownloader(date=date, revision=revision)\n",
    "    future = client.submit(\n",
    "        func=source.discover, key=f\"download-{date}\"\n",
    "    )  # triggers download of the file(s), or loads from cache\n",
    "    futures.append(future)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 64/64 [2:24:35<00:00, 135.56s/it]  \n"
     ]
    }
   ],
   "source": [
    "# Check download progress here, https://stackoverflow.com/a/37901797/6611055\n",
    "responses = []\n",
    "for f in tqdm.tqdm(\n",
    "    iterable=dask.distributed.as_completed(futures=futures), total=len(futures)\n",
    "):\n",
    "    responses.append(f.result())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In case of error, check which downloads are unfinished\n",
    "# Manually delete those folders and retry\n",
    "unfinished = []\n",
    "for foo in futures:\n",
    "    if foo.status != \"finished\":\n",
    "        print(foo)\n",
    "        unfinished.append(foo)\n",
    "        # foo.retry()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    assert len(unfinished) == 0\n",
    "except AssertionError:\n",
    "    for task in unfinished:\n",
    "        print(task)\n",
    "    raise ValueError(\n",
    "        f\"{len(unfinished)} download tasks are unfinished,\"\n",
    "        \" please delete those folders and retry again!\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploratory data analysis on local files\n",
    "\n",
    "Now that we've downloaded a good chunk of data and cached them locally,\n",
    "we can have some fun with visualizing the point clouds!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "root_directory = os.path.dirname(\n",
    "    catalog.icesat2atl06.storage_options[\"simplecache\"][\"cache_storage\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "def get_crossing_dates(\n",
    "    catalog_entry: intake.catalog.local.LocalCatalogEntry,\n",
    "    root_directory: str,\n",
    "    referencegroundtrack: str = \"????\",\n",
    "    datetimestr: str = \"*\",\n",
    "    cyclenumber: str = \"??\",\n",
    "    orbitalsegment: str = \"??\",\n",
    "    version: str = \"003\",\n",
    "    revision: str = \"01\",\n",
    ") -> dict:\n",
    "    \"\"\"\n",
    "    Given a 4-digit reference groundtrack (e.g. 1234),\n",
    "    we output a dictionary where the\n",
    "    key is the date in \"YYYY.MM.DD\" format when an ICESAT2 crossing was made and the\n",
    "    value is the filepath to the HDF5 data file.\n",
    "    \"\"\"\n",
    "\n",
    "    # Get a glob string that looks like \"ATL06_??????????????_XXXX????_002_01.h5\"\n",
    "    globpath: str = catalog_entry.path_as_pattern\n",
    "    if datetimestr == \"*\":\n",
    "        globpath: str = globpath.replace(\"{datetime:%Y%m%d%H%M%S}\", \"??????????????\")\n",
    "    globpath: str = globpath.format(\n",
    "        referencegroundtrack=referencegroundtrack,\n",
    "        cyclenumber=cyclenumber,\n",
    "        orbitalsegment=orbitalsegment,\n",
    "        version=version,\n",
    "        revision=revision,\n",
    "    )\n",
    "\n",
    "    # Get list of filepaths (dates are contained in the filepath)\n",
    "    globedpaths: list = glob.glob(os.path.join(root_directory, \"??????????\", globpath))\n",
    "\n",
    "    # Pick out just the dates in \"YYYY.MM.DD\" format from the globedpaths\n",
    "    # crossingdates = [os.path.basename(os.path.dirname(p=p)) for p in globedpaths]\n",
    "    crossingdates: dict = {\n",
    "        os.path.basename(os.path.dirname(p=p)): p for p in sorted(globedpaths)\n",
    "    }\n",
    "\n",
    "    return crossingdates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "crossing_dates_dict = {}\n",
    "for rgt in range(1, 1388):  # ReferenceGroundTrack goes from 0001 to 1387\n",
    "    referencegroundtrack: str = f\"{rgt}\".zfill(4)\n",
    "    crossing_dates: dict = dask.delayed(get_crossing_dates)(\n",
    "        catalog_entry=catalog.icesat2atl06,\n",
    "        root_directory=root_directory,\n",
    "        referencegroundtrack=referencegroundtrack,\n",
    "    )\n",
    "    crossing_dates_dict[referencegroundtrack] = crossing_dates\n",
    "crossing_dates_dict = dask.compute(crossing_dates_dict)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['2018.10.21', '2019.01.20', '2019.04.21', '2019.10.19', '2020.01.18', '2020.04.18'])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "crossing_dates_dict[\"0349\"].keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![ICESat-2 Laser Beam Pattern](https://ars.els-cdn.com/content/image/1-s2.0-S0034425719303712-gr1.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def six_laser_beams(filepaths: list) -> dask.dataframe.DataFrame:\n",
    "    \"\"\"\n",
    "    For all 6 lasers along one reference ground track,\n",
    "    concatenate all points from all crossing dates into one Dask DataFrame\n",
    "\n",
    "    E.g. if there are 5 crossing dates and 6 lasers,\n",
    "    there would be data from 5 x 6 = 30 files being concatenated together.\n",
    "    \"\"\"\n",
    "    lasers: list = [\"gt1l\", \"gt1r\", \"gt2l\", \"gt2r\", \"gt3l\", \"gt3r\"]\n",
    "\n",
    "    objs: list = [\n",
    "        xr.open_mfdataset(\n",
    "            paths=filepaths,\n",
    "            combine=\"by_coords\",\n",
    "            engine=\"h5netcdf\",\n",
    "            group=f\"{laser}/land_ice_segments\",\n",
    "            parallel=True,\n",
    "        ).assign_coords(coords={\"laser\": laser})\n",
    "        for laser in lasers\n",
    "    ]\n",
    "\n",
    "    try:\n",
    "        da: xr.Dataset = xr.concat(objs=objs, dim=\"laser\")\n",
    "        df: dask.dataframe.DataFrame = da.unify_chunks().to_dask_dataframe()\n",
    "    except ValueError:\n",
    "        # ValueError: cannot reindex or align along dimension 'delta_time'\n",
    "        # because the index has duplicate values\n",
    "        df: dask.dataframe.DataFrame = dask.dataframe.concat(\n",
    "            [obj.unify_chunks().to_dask_dataframe() for obj in objs]\n",
    "        )\n",
    "\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset_dict = {}\n",
    "# ReferenceGroundTrack goes from 0001 to 1387\n",
    "for referencegroundtrack in list(crossing_dates_dict)[348:349]:\n",
    "    # print(referencegroundtrack)\n",
    "    filepaths = list(crossing_dates_dict[referencegroundtrack].values())\n",
    "    if len(filepaths) > 0:\n",
    "        dataset_dict[referencegroundtrack] = dask.delayed(obj=six_laser_beams)(\n",
    "            filepaths=filepaths\n",
    "        )\n",
    "        # df = six_laser_beams(filepaths=filepaths)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = dataset_dict[\"0349\"].compute()  # loads into a dask dataframe (lazy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><strong>Dask DataFrame Structure:</strong></div>\n",
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>delta_time</th>\n",
       "      <th>laser</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>atl06_quality_summary</th>\n",
       "      <th>h_li</th>\n",
       "      <th>h_li_sigma</th>\n",
       "      <th>segment_id</th>\n",
       "      <th>sigma_geo_h</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>npartitions=190</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>datetime64[ns]</td>\n",
       "      <td>object</td>\n",
       "      <td>float64</td>\n",
       "      <td>float64</td>\n",
       "      <td>float64</td>\n",
       "      <td>float32</td>\n",
       "      <td>float32</td>\n",
       "      <td>float64</td>\n",
       "      <td>float32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>340326</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10697850</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11087807</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>\n",
       "<div>Dask Name: concat-indexed, 26334 tasks</div>"
      ],
      "text/plain": [
       "Dask DataFrame Structure:\n",
       "                     delta_time   laser latitude longitude atl06_quality_summary     h_li h_li_sigma segment_id sigma_geo_h\n",
       "npartitions=190                                                                                                            \n",
       "0                datetime64[ns]  object  float64   float64               float64  float32    float32    float64     float32\n",
       "340326                      ...     ...      ...       ...                   ...      ...        ...        ...         ...\n",
       "...                         ...     ...      ...       ...                   ...      ...        ...        ...         ...\n",
       "10697850                    ...     ...      ...       ...                   ...      ...        ...        ...         ...\n",
       "11087807                    ...     ...      ...       ...                   ...      ...        ...        ...         ...\n",
       "Dask Name: concat-indexed, 26334 tasks"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute every referencegroundtrack, slow... though somewhat parallelized\n",
    "# dataset_dict = dask.compute(dataset_dict)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# big dataframe containing data across all 1387 reference ground tracks!\n",
    "# bdf = dask.dataframe.concat(dfs=list(dataset_dict.values()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# https://xarray.pydata.org/en/stable/combining.html#concatenate\n",
    "# For all 6 lasers one one date ~~along one reference ground track~~,\n",
    "# concatenate all points ~~from one dates~~ into one xr.Dataset\n",
    "lasers = [\"gt1l\", \"gt1r\", \"gt2l\", \"gt2r\", \"gt3l\", \"gt3r\"]\n",
    "da = xr.concat(\n",
    "    objs=(\n",
    "        catalog.icesat2atl06(laser=laser, referencegroundtrack=referencegroundtrack)\n",
    "        .to_dask()\n",
    "        for laser in lasers\n",
    "    ),\n",
    "    dim=pd.Index(data=lasers, name=\"laser\")\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot ATL06 points!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert dask.DataFrame to pd.DataFrame\n",
    "df: pd.DataFrame = df.compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop points with poor quality\n",
    "df = df.dropna(subset=[\"h_li\"]).query(expr=\"atl06_quality_summary == 0\").reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>index</th>\n",
       "      <th>delta_time</th>\n",
       "      <th>laser</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>atl06_quality_summary</th>\n",
       "      <th>h_li</th>\n",
       "      <th>h_li_sigma</th>\n",
       "      <th>segment_id</th>\n",
       "      <th>sigma_geo_h</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>800772</th>\n",
       "      <td>4831606</td>\n",
       "      <td>2019-04-21 03:46:39.425294928</td>\n",
       "      <td>gt3l</td>\n",
       "      <td>-76.712393</td>\n",
       "      <td>50.365090</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3448.515625</td>\n",
       "      <td>0.009265</td>\n",
       "      <td>1578380.0</td>\n",
       "      <td>0.307556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1524633</th>\n",
       "      <td>9459271</td>\n",
       "      <td>2020-04-18 10:26:05.075350352</td>\n",
       "      <td>gt1r</td>\n",
       "      <td>-75.535912</td>\n",
       "      <td>49.316946</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3399.856934</td>\n",
       "      <td>0.010297</td>\n",
       "      <td>1585119.0</td>\n",
       "      <td>0.345309</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1241327</th>\n",
       "      <td>7669491</td>\n",
       "      <td>2020-01-18 14:46:17.080013576</td>\n",
       "      <td>gt2r</td>\n",
       "      <td>-75.455553</td>\n",
       "      <td>49.387659</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3389.260742</td>\n",
       "      <td>0.010744</td>\n",
       "      <td>1585549.0</td>\n",
       "      <td>0.319406</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1386790</th>\n",
       "      <td>8632127</td>\n",
       "      <td>2020-01-18 14:48:24.751603336</td>\n",
       "      <td>gt3r</td>\n",
       "      <td>-67.439814</td>\n",
       "      <td>46.033186</td>\n",
       "      <td>0.0</td>\n",
       "      <td>26.925806</td>\n",
       "      <td>0.013543</td>\n",
       "      <td>1630850.0</td>\n",
       "      <td>0.301213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1681485</th>\n",
       "      <td>10423283</td>\n",
       "      <td>2020-04-18 10:27:26.745665840</td>\n",
       "      <td>gt3r</td>\n",
       "      <td>-70.406813</td>\n",
       "      <td>47.045861</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2529.187256</td>\n",
       "      <td>0.019172</td>\n",
       "      <td>1614106.0</td>\n",
       "      <td>0.301355</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            index                    delta_time laser   latitude  longitude  \\\n",
       "800772    4831606 2019-04-21 03:46:39.425294928  gt3l -76.712393  50.365090   \n",
       "1524633   9459271 2020-04-18 10:26:05.075350352  gt1r -75.535912  49.316946   \n",
       "1241327   7669491 2020-01-18 14:46:17.080013576  gt2r -75.455553  49.387659   \n",
       "1386790   8632127 2020-01-18 14:48:24.751603336  gt3r -67.439814  46.033186   \n",
       "1681485  10423283 2020-04-18 10:27:26.745665840  gt3r -70.406813  47.045861   \n",
       "\n",
       "         atl06_quality_summary         h_li  h_li_sigma  segment_id  \\\n",
       "800772                     0.0  3448.515625    0.009265   1578380.0   \n",
       "1524633                    0.0  3399.856934    0.010297   1585119.0   \n",
       "1241327                    0.0  3389.260742    0.010744   1585549.0   \n",
       "1386790                    0.0    26.925806    0.013543   1630850.0   \n",
       "1681485                    0.0  2529.187256    0.019172   1614106.0   \n",
       "\n",
       "         sigma_geo_h  \n",
       "800772      0.307556  \n",
       "1524633     0.345309  \n",
       "1241327     0.319406  \n",
       "1386790     0.301213  \n",
       "1681485     0.301355  "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get a small random sample of our data\n",
    "dfs = df.sample(n=1_000, random_state=42)\n",
    "dfs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs.hvplot.scatter(\n",
    "    x=\"longitude\",\n",
    "    y=\"latitude\",\n",
    "    by=\"laser\",\n",
    "    hover_cols=[\"delta_time\", \"segment_id\"],\n",
    "    # datashade=True, dynspread=True,\n",
    "    # width=800, height=500, colorbar=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transform from EPSG:4326 (lat/lon) to EPSG:3031 (Antarctic Polar Stereographic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs[\"x\"], dfs[\"y\"] = deepicedrain.lonlat_to_xy(\n",
    "    longitude=dfs.longitude, latitude=dfs.latitude\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>index</th>\n",
       "      <th>delta_time</th>\n",
       "      <th>laser</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>atl06_quality_summary</th>\n",
       "      <th>h_li</th>\n",
       "      <th>h_li_sigma</th>\n",
       "      <th>segment_id</th>\n",
       "      <th>sigma_geo_h</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>800772</th>\n",
       "      <td>4831606</td>\n",
       "      <td>2019-04-21 03:46:39.425294928</td>\n",
       "      <td>gt3l</td>\n",
       "      <td>-76.712393</td>\n",
       "      <td>50.365090</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3448.515625</td>\n",
       "      <td>0.009265</td>\n",
       "      <td>1578380.0</td>\n",
       "      <td>0.307556</td>\n",
       "      <td>1.116662e+06</td>\n",
       "      <td>9.249294e+05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1524633</th>\n",
       "      <td>9459271</td>\n",
       "      <td>2020-04-18 10:26:05.075350352</td>\n",
       "      <td>gt1r</td>\n",
       "      <td>-75.535912</td>\n",
       "      <td>49.316946</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3399.856934</td>\n",
       "      <td>0.010297</td>\n",
       "      <td>1585119.0</td>\n",
       "      <td>0.345309</td>\n",
       "      <td>1.197869e+06</td>\n",
       "      <td>1.029714e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1241327</th>\n",
       "      <td>7669491</td>\n",
       "      <td>2020-01-18 14:46:17.080013576</td>\n",
       "      <td>gt2r</td>\n",
       "      <td>-75.455553</td>\n",
       "      <td>49.387659</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3389.260742</td>\n",
       "      <td>0.010744</td>\n",
       "      <td>1585549.0</td>\n",
       "      <td>0.319406</td>\n",
       "      <td>1.205870e+06</td>\n",
       "      <td>1.034007e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1386790</th>\n",
       "      <td>8632127</td>\n",
       "      <td>2020-01-18 14:48:24.751603336</td>\n",
       "      <td>gt3r</td>\n",
       "      <td>-67.439814</td>\n",
       "      <td>46.033186</td>\n",
       "      <td>0.0</td>\n",
       "      <td>26.925806</td>\n",
       "      <td>0.013543</td>\n",
       "      <td>1630850.0</td>\n",
       "      <td>0.301213</td>\n",
       "      <td>1.786477e+06</td>\n",
       "      <td>1.723183e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1681485</th>\n",
       "      <td>10423283</td>\n",
       "      <td>2020-04-18 10:27:26.745665840</td>\n",
       "      <td>gt3r</td>\n",
       "      <td>-70.406813</td>\n",
       "      <td>47.045861</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2529.187256</td>\n",
       "      <td>0.019172</td>\n",
       "      <td>1614106.0</td>\n",
       "      <td>0.301355</td>\n",
       "      <td>1.572852e+06</td>\n",
       "      <td>1.464356e+06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            index                    delta_time laser   latitude  longitude  \\\n",
       "800772    4831606 2019-04-21 03:46:39.425294928  gt3l -76.712393  50.365090   \n",
       "1524633   9459271 2020-04-18 10:26:05.075350352  gt1r -75.535912  49.316946   \n",
       "1241327   7669491 2020-01-18 14:46:17.080013576  gt2r -75.455553  49.387659   \n",
       "1386790   8632127 2020-01-18 14:48:24.751603336  gt3r -67.439814  46.033186   \n",
       "1681485  10423283 2020-04-18 10:27:26.745665840  gt3r -70.406813  47.045861   \n",
       "\n",
       "         atl06_quality_summary         h_li  h_li_sigma  segment_id  \\\n",
       "800772                     0.0  3448.515625    0.009265   1578380.0   \n",
       "1524633                    0.0  3399.856934    0.010297   1585119.0   \n",
       "1241327                    0.0  3389.260742    0.010744   1585549.0   \n",
       "1386790                    0.0    26.925806    0.013543   1630850.0   \n",
       "1681485                    0.0  2529.187256    0.019172   1614106.0   \n",
       "\n",
       "         sigma_geo_h             x             y  \n",
       "800772      0.307556  1.116662e+06  9.249294e+05  \n",
       "1524633     0.345309  1.197869e+06  1.029714e+06  \n",
       "1241327     0.319406  1.205870e+06  1.034007e+06  \n",
       "1386790     0.301213  1.786477e+06  1.723183e+06  \n",
       "1681485     0.301355  1.572852e+06  1.464356e+06  "
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dfs.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dfs.hvplot.scatter(\n",
    "    x=\"x\",\n",
    "    y=\"y\",\n",
    "    by=\"laser\",\n",
    "    hover_cols=[\"delta_time\", \"segment_id\", \"h_li\"],\n",
    "    # datashade=True, dynspread=True,\n",
    "    # width=800, height=500, colorbar=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot cross section view\n",
    "dfs.hvplot.scatter(x=\"x\", y=\"h_li\", by=\"laser\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experimental Work-in-Progress stuff below"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Old way of making a DEM grid surface from points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://github.com/ICESAT-2HackWeek/gridding/blob/master/notebook/utils.py#L23\n",
    "def make_grid(xmin, xmax, ymin, ymax, dx, dy):\n",
    "    \"\"\"Construct output grid-coordinates.\"\"\"\n",
    "\n",
    "    # Setup grid dimensions\n",
    "    Nn = int((np.abs(ymax - ymin)) / dy) + 1\n",
    "    Ne = int((np.abs(xmax - xmin)) / dx) + 1\n",
    "\n",
    "    # Initiate x/y vectors for grid\n",
    "    x_i = np.linspace(xmin, xmax, num=Ne)\n",
    "    y_i = np.linspace(ymin, ymax, num=Nn)\n",
    "\n",
    "    return np.meshgrid(x_i, y_i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xi, yi = make_grid(\n",
    "    xmin=dfs.x.min(), xmax=dfs.x.max(), ymin=dfs.y.max(), ymax=dfs.y.min(), dx=10, dy=10\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ar = scipy.interpolate.griddata(points=(dfs.x, dfs.y), values=dfs.h_li, xi=(xi, yi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(ar, extent=(dfs.x.min(), dfs.x.max(), dfs.y.min(), dfs.y.max()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import plotly.express as px"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "px.scatter_3d(data_frame=dfs, x=\"longitude\", y=\"latitude\", z=\"h_li\", color=\"laser\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Play using XrViz\n",
    "\n",
    "Install the PyViz JupyterLab extension first using the [extension manager](https://jupyterlab.readthedocs.io/en/stable/user/extensions.html#using-the-extension-manager) or via the command below:\n",
    "\n",
    "```bash\n",
    "jupyter labextension install @pyviz/jupyterlab_pyviz@v0.8.0 --no-build\n",
    "jupyter labextension list  # check to see that extension is installed\n",
    "jupyter lab build --debug  # build extension ??? with debug messages printed\n",
    "```\n",
    "\n",
    "Note: Had to add `network-timeout 600000` to `.yarnrc` file to resolve university network issues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xrviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xrviz.example()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://xrviz.readthedocs.io/en/latest/set_initial_parameters.html\n",
    "initial_params = {\n",
    "    # Select variable to plot\n",
    "    \"Variables\": \"h_li\",\n",
    "    # Set coordinates\n",
    "    \"Set Coords\": [\"longitude\", \"latitude\"],\n",
    "    # Axes\n",
    "    \"x\": \"longitude\",\n",
    "    \"y\": \"latitude\",\n",
    "    # \"sigma\": \"animate\",\n",
    "    # Projection\n",
    "    # \"is_geo\": True,\n",
    "    # \"basemap\": True,\n",
    "    # \"crs\": \"PlateCarree\"\n",
    "}\n",
    "dashboard = xrviz.dashboard.Dashboard(data=dataset)  # , initial_params=initial_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dashboard.panel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dashboard.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OpenAltimetry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"minx=-154.56678505984297&miny=-88.82881451427136&maxx=-125.17872921546498&maxy=-81.34051361301398&date=2019-05-02&trackId=516\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paste the OpenAltimetry selection parameters here\n",
    "OA_REFERENCE_URL = \"minx=-177.64275595145213&miny=-88.12014866942751&maxx=-128.25920892322736&maxy=-85.52394234080862&date=2019-05-02&trackId=515\"\n",
    "# We populate a list with the photon data using the OpenAltimetry API, no HDF!\n",
    "OA_URL = (\n",
    "    \"https://openaltimetry.org/data/icesat2/getPhotonData?client=jupyter&\"\n",
    "    + OA_REFERENCE_URL\n",
    ")\n",
    "OA_PHOTONS = [\"Noise\", \"Low\", \"Medium\", \"High\"]\n",
    "# OA_PLOTTED_BEAMS = [1,2,3,4,5,6] you can select up to 6 beams for each ground track.\n",
    "# Some beams may not be usable due cloud covering or QC issues.\n",
    "OA_BEAMS = [3, 4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minx, miny, maxx, maxy = [-156, -88, -127, -84]\n",
    "date = \"2019-05-02\"  # UTC date?\n",
    "track = 515  #\n",
    "beam = 1  # 1 to 6\n",
    "params = {\n",
    "    \"client\": \"jupyter\",\n",
    "    \"minx\": minx,\n",
    "    \"miny\": miny,\n",
    "    \"maxx\": maxx,\n",
    "    \"maxy\": maxy,\n",
    "    \"date\": date,\n",
    "    \"trackId\": str(track),\n",
    "    \"beam\": str(beam),\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = requests.get(\n",
    "    url=\"https://openaltimetry.org/data/icesat2/getPhotonData\", params=params\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OpenAltimetry Data cleansing\n",
    "df = pd.io.json.json_normalize(data=r.json()[\"series\"], meta=\"name\", record_path=\"data\")\n",
    "df.name = df.name.str.split().str.get(0)  # Get e.g. just \"Low\" instead of \"Low [12345]\"\n",
    "df.query(\n",
    "    expr=\"name in ('Low', 'Medium', 'High')\", inplace=True\n",
    ")  # filter out Noise and Buffer points\n",
    "\n",
    "df.rename(columns={0: \"latitude\", 1: \"elevation\", 2: \"longitude\"}, inplace=True)\n",
    "df = df.reindex(\n",
    "    columns=[\"longitude\", \"latitude\", \"elevation\", \"name\"]\n",
    ")  # reorder columns\n",
    "df.reset_index(inplace=True)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.hvplot.scatter(x=\"latitude\", y=\"elevation\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,py:hydrogen"
  },
  "kernelspec": {
   "display_name": "deepicedrain",
   "language": "python",
   "name": "deepicedrain"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
